不一定正确的多分组差异分析结果热图展现

最近给学徒安排了一个文献数据挖掘任务，但是检查她结果发现跟文献差异有点大，因为时间关系，来不及全盘细致检查，所以先发出来交给读者考量。主要是基于我前些天在生信技能树分享的一个策略：如果你的分组比较多，差异分析策略有哪些？  

1.复现的图片

Paper:Machine learning analysis of gene expression data reveals novel diagnostic and prognostic biomarkers and identifies therapeutic targets for soft tissue sarcomas

2.术语：

3.Something to clarify

• 不是直接从TCGA下载的，而是文章整理的，文章是下面这篇：

Comprehensive and Integrated Genomic Characterization of Adult Soft Tissue Sarcomas

• 这里的热图不是用表达矩阵作图，而是，用的logFC；

• 低表达基因作过滤用的是，cpm<2，差异分析用的是logFC>0和Adjust.pvalue <0.05；

• GO富集分析用的是enrichR；

Step01-Data download

• GDCquery+GDCdownload+GDCprepare进行数据下载, 大家可能需要从头理解TCGA背景知识：悄咪咪的上线了TCGA知识图谱视频教程（B站和YouTube直达）

rm(list=ls())
library(TCGAbiolinks)
library(SummarizedExperiment)
query <- GDCquery(project = "TCGA-SARC",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")
GDCdownload(query)
SARC_DATA<-GDCprepare(query,save=FALSE)
save(SARC_DATA,file='sarc.Rdata')

Step02-Data_for_analyse

• 上文有写，临床信息来自一篇文章，NIHMS912614-supplement-10.xlsx为其附表；

rm(list=ls())
load('sarc.Rdata')
library(SummarizedExperiment)
summary(SARC_DATA)
sarc_data<- assay(SARC_DATA)
dim(sarc_data)
####临床信息来源于另一篇paper的附表
library('openxlsx')
clinic <- read.xlsx('NIHMS912614-supplement-10.xlsx',sheet = 2,startRow = 2)
colnames(sarc_data)<- substr(colnames(sarc_data),start = 1,stop = 15)
sarc_choose<- sarc_data[,clinic$TCGA.barcode]
save(sarc_choose,clinic,file='ex_pd.Rdata')

Step03-Differential analysis

• 阈值设定（cpm）均源于我对文章的理解，差异分析的步骤也是如此；

• 这里整个差异分析过程，采用了循环的方式进行，里面主要在循环的就是group_list和design，可能需要详细看才能理解；

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'ex_pd.Rdata')
library(edgeR)
#####过滤掉cpm<2的counts
tmp<- cpm(sarc_choose)
tmp1<- tmp[rowSums(tmp<2)>1,]
sarc_choose<- sarc_choose[rownames(tmp1),]
identical(colnames(sarc_choose),clinic$TCGA.barcode)
####用TMM方法进行normalization
dge <- DGEList(counts=sarc_choose)
dge <- calcNormFactors(dge,method = 'TMM')
####构建对应的group_list
subtype <- unique(clinic$short.histo)
group_list <- c()
for( i in 1:length(subtype)){
  tmp<- ifelse(clinic$short.histo==subtype[i],subtype[i],'other')
  group_list<- cbind(group_list,tmp)
  colnames(group_list)[i] <- subtype[i]
}
#####差异分析和PCA 
deg_re <- matrix()
for( i in 1:ncol(group_list)){
  design <- model.matrix(~0+factor(group_list[,i]))
  colnames(design)=levels(factor(group_list[,i]))
  rownames(design)=colnames(sarc_choose)
  data_v <- voom(dge,design)
  ###PCA
  library("FactoMineR")
  library("factoextra") 
  ####
  df.pca <- PCA(t(data_v$E), graph = FALSE)
  fviz_pca_ind(df.pca,
               geom.ind = "point",
               col.ind = group_list[,i], 
               addEllipses = TRUE, 
               legend.title = "Groups"
  )
  title<- paste0('./figures/',colnames(group_list)[i],'.png')
  ggsave(title)
  constrsts<- paste0(colnames(group_list)[i],"-other")
  contrast.matrix<-makeContrasts(contrasts=constrsts,levels = design)
  ######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
  fit <- lmFit(data_v,design)
  fit1 <- contrasts.fit(fit, contrast.matrix)
  fit2 <- eBayes(fit1, 0.01)
  ####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
  tT_tmp <- topTable(fit2, adjust="fdr", number=nrow(fit2))
  tT_tmp<- subset(tT_tmp, select=c('adj.P.Val',"P.Value","logFC"))
  colnames(tT_tmp)<- paste0(colnames(group_list)[i],colnames(tT_tmp))
  deg_re <- cbind(deg_re,tT_tmp)
  }
deg_re <- deg_re[,-1]
save(deg_re,file='deg.Rdata')

Step04-pheatmap

• 文章是用logFC做热图，这里等于挑选了，差异基因并集的logFC；

• 之前直接用logFC>0的阈值进行热图绘制，热图很丑，后面把阈值调高了进行绘图（曾老师原话是排序，这里并没有这样）；

• 看文章的图，像是聚类后的，但又没有展示tree，我就把聚类结果存为变量，再进行的热图绘制；

rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
histo<- unique(clinic$short.histo)
keep <- c()
for(i in 1:length(histo)){
  tT<- deg_re[,grepl(histo[i],colnames(deg_re))]
  colnames(tT)<- gsub(histo[i],'',colnames(tT))
  ####这个logFC的阈值是文章选的
  tT$change = ifelse(tT$adj.P.Val < 0.05 & tT$logFC>2,TRUE,FALSE)
  keep <- cbind(keep,tT$change)
}
#####文章是用logFC做热图，这里等于挑选了，差异基因并集的logFC
logFC_pheatmap<- deg_re[apply(keep,1,function(x){sum(x)>0}),grepl('logFC',colnames(deg_re))]
colnames(logFC_pheatmap) <- gsub('logFC','',colnames(logFC_pheatmap))
nrow(logFC_pheatmap)
library(pheatmap)
tmp<- pheatmap(logFC_pheatmap,silent=T)
####调出热图聚类结果，并存储变量
final_phe<- logFC_pheatmap[tmp$tree_row$order,tmp$tree_col$order]
bk = c(seq(min(final_phe),max(final_phe), length=100))
pheatmap(final_phe,
         breaks=bk,
         show_rownames = F,
         cluster_rows=F,
         cluster_cols = F,
         color = c(colorRampPalette(c("blue", 'white','red'))(100)))

Step05-GO enrichment

• paper中富集所用的database为GO_Biological_Process_2015和R包为enrichR；

• 最后的富集结果，某些相应的Adjust.pvalue 没有满足条件的，我选了前三条，共7+7*3行，这里截图展示；

rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
histo<- unique(clinic$short.histo)
library(org.Hs.eg.db)
E2S<-select(org.Hs.eg.db,keys=keys(org.Hs.eg.db),columns = c("SYMBOL",'ENSEMBL'))
deg_re$ENSEMBL <- rownames(deg_re)
deg_re<- merge(deg_re,E2S,by='ENSEMBL')
deg_re <- na.omit(deg_re)
GO_histo <- list()
for(i in 1:length(histo)){
  tT<- deg_re[,grepl(histo[i],colnames(deg_re))]
  colnames(tT)<- gsub(histo[i],'',colnames(tT))
  tT$change = ifelse(tT$adj.P.Val < 0.05 & tT$logFC0,TRUE,FALSE)
######paper中富集所用的database为GO_Biological_Process_2015和R包为enrichR
  library(enrichR)
  dbs <- c("GO_Biological_Process_2015")
  enriched <- enrichr(deg_re$SYMBOL[tT$change], dbs)
  bp <- enriched[["GO_Biological_Process_2015"]]
  GO_histo[[histo[i]]]<- bp[1:3,]
}
GO_final <- c()
for(i in 1:length(GO_histo)){
  name <- paste0(names(GO_histo)[i])
  tmp<- GO_histo[[i]]
  tmp<- c(name,as.character(tmp[,grepl('Term',colnames(tmp))]))
  GO_final <- c(GO_final,tmp)
}
GO_final <- as.data.frame(GO_final)

写在后面

因为并没有看到GO的数据库注释结果，而且在我检查热图里面各个癌症亚型的差异基因数量的时候，发现跟原文不一致！

首先，可以看到原文数量如下:

然后我使用代码检查学徒的差异分析结果：

rm(list=ls())
load('deg.Rdata')
load('ex_pd.Rdata')
sarc_choose=sarc_choose[rownames(sarc_choose) %in% rownames(deg_re),]
histo<- unique(clinic$short.histo)

keep=do.call(cbind,lapply(histo, function(x){
  print(x)
  cc=deg_re[,grep(x,colnames(deg_re))]
  return(cc[,1]<0.05 & cc[,3]>2)
}))
colnames(keep)=histo
rownames(keep)=rownames(deg_re)
apply(keep,2,table)

得到的差异基因数量分布如下：

原文是DDLPS的差异基因数量最多，而学徒代码里面是 SS 的最多，所以中间过程还是需要仔细检查，争取找到差异。

不过，今天是国庆节最后一天，我还在路上，就不耗费时间在这上面了，感兴趣的粉丝可以根据推文里面的代码来尝试学习+理解+发给我你的见解！

最后的最后

我还想说一下的是，其实这个项目的样本数量已经超过200，而且分组是7个，已经可以使用单细胞数据分析策略了，比如我这里就使用seurat对表达矩阵进行tSNE可视化：

还可以使用分组，看看与文章的相似性

文章分组如下：

还有更多单细胞探索，大家自行学习，比如 DoHeatmap 可以绘制我给学徒的任务。

号外：（南京、南宁见！）全国巡讲第17-18站（生信入门课加量不加价）